Location determination using doppler and pseudorange measurements from fewer than four satellites

ABSTRACT

A method for estimating two or three location coordinates and/or receiver time offset for a receiver that receives and analyzes timed signals transmitted from one, two or three non-geosynchronous satellites (GPS, GLONASS, modified LEO, etc.). Combinations of pseudorange, Doppler shift and/or first time derivative and/or second time derivative of Doppler shift are determined and used to estimate the location coordinates and/or receiver time offset.

FIELD OF THE INVENTION

This invention relates to use of signals received from fewer than four satellites to determine location of a signal receiver.

BACKGROUND OF THE INVENTION

Location determination (LD) systems, such as the Global Positioning System (GPS), the Global Orbiting Navigational Satellite System (GLONASS) and the Low Earth Orbit (LEO) System, receive and analyze timed signals from a plurality of satellites (usually three or more) in order to estimate the location and time offset of an LD signal receiver. The LD signals received may be used to measure pseudoranges from the receiver to each satellite, carrier phases, signal Doppler shifts and the like, but it is usually necessary to receive and analyze these signals from many satellites in order to estimate the location coordinates and time offset or related quantities.

In some situations, acceptable signals from an adequate number of satellites is not available, because of signal interference from an adjacent urban canyon or other structure(s), because of an unfortunate geometric configuration for the contributing satellites, because of partial or full masking of signals transmitted near the horizon, or because of non-availability of signal processing in one or more receiver channels. What is needed is a signal receiver system that can receive and analyze LD signals from fewer than four satellites to estimate location coordinates and receiver time offset, or other relevant variables. Preferably, the system should work with LD signals received from as few as one or two LD signal sources (usually, but not necessarily, satellites) at substantially the same signal receipt time. Preferably, the system should work with an arbitrary satellite configuration, not merely with satellites in equatorial orbits.

SUMMARY OF THE INVENTION

These needs are met by the invention, which provides method and apparatus that receives LD signals from fewer than four LD signal sources, measures pseudorange, signal Doppler shift and selected time derivatives of Doppler shift, and estimates the location coordinates and time offset ("location fix coordinates") for the LD signal receiver. In one approach, pseudorange and Doppler shift (no time derivatives) for each of two LD signal sources are measured and analyzed to estimate the location fix coordinates. In another approach, pseudorange and Doppler shift (no time derivatives) for a single moving LD signal source are measured and analyzed at two distinct times to estimate the location fix coordinates. In another approach, any two of the pseudorange and the zeroth, first and second time derivatives of Doppler shift at two distinct times are measured and analyzed to estimate the location fix coordinates. In another approach, pseudorange and at least one of the zeroth, first and second time derivatives of Doppler shift for a single moving LD signal source are measured and analyzed at two distinct times to estimate the location fix coordinates. In another approach, pseudorange, Doppler shift and the first and second time derivatives of Doppler shift are measured and analyzed at a selected time to estimate the location fix coordinates. In another approach, any combination of four of the pseudorange and the zeroth, first and second time derivatives of Doppler shift at two or more distinct times are measured and analyzed to estimate the location fix coordinates. If one or more of the location fix coordinates (x_(r),y_(r),z_(r),b_(r)) is known, the number of measurements is reduced accordingly.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 illustrates a situation in which the invention can be used.

FIGS. 2-4 are flow charts illustrating methods to practice the invention.

FIG. 5 illustrates analysis of an elliptical satellite orbit.

FIG. 6 illustrates apparatus to practice the invention.

DESCRIPTION OF BEST MODES OF THE INVENTION

In FIG. 1, an LD signal antenna 12 and receiver 11 (referred to collectively as an "LD receiver") receives LD signals from one or more moving LD signal sources 13A, 13B, 13C, whose location coordinates are known as a function of time t and which move relative to the receiver. The LD receiver 11 rotates with the Earth's surface, has initial location coordinates (x_(r0),y_(r0),z_(r0)), and has fixed receiver time offset Δt_(r) with b_(r) =c·Δt_(r), where c is the velocity of light in the ambient medium. The location fix coordinates (x_(r),y_(r),z_(r),b_(r)) for the LD receiver 11 are to be estimated, using LD signals received from one, two or three LD signal sources.

The LD signal receiver 11 occupies a fixed location relative to the Earth's surface as the Earth rotates on its axis. Expressed as a function of time relative to an absolute coordinate system in space, the location coordinates for the LD receiver 11 can be expressed in a Cartesian system, with origin at the Earth's center of rotation and z-axis coinciding with the Earth's rotation axis, as

    x.sub.r (t)=x.sub.r0 cos ω.sub.e (t-t0)+y.sub.r0 sin ω.sub.e (t-t0),                                                   (1)

    y.sub.r (t)=-x.sub.r0 sin ω.sub.e (t-t0)+y.sub.r0 cos ω.sub.e (t-t0),                                                   (2)

    z.sub.r (t)=z.sub.r0,                                      (3)

where x_(r0), y_(r0) and z_(r0) are the (unknown) location coordinate values of the LD receiver at a selected time t'=t-t0 and ω_(e) is the known angular velocity of the Earth.

The first, second and third time derivatives of the location coordinates for the LD signal receiver 11 are found to be

    ∂x.sub.r /∂t=ω.sub.e {-x.sub.r0 sin ω.sub.e (t-t0)+y.sub.r0 cos ω.sub.e (t-t0)},  (4)

    ∂y.sub.r /∂t=ω.sub.e {-x.sub.r0 cos ω.sub.e (t-t0)-y.sub.r0 sin ω.sub.e (t-t0)},  (5)

    ∂z.sub.r /∂t=0,                  (6)

    ∂.sup.2 x.sub.r /∂t.sup.2 =ω.sub.e.sup.2 {-x.sub.r0 cos ω.sub.e (t-t0)-y.sub.r0 sin ω.sub.e (t-t0)},(7)

    ∂.sup.2 x.sub.r /∂t.sup.2 =ω.sub.e.sup.2 {-x.sub.r0 sin ω.sub.e (t-t0)-y.sub.r0 cos ω.sub.e (t-t0)},(8)

    ∂.sup.2 z.sub.r /∂t.sup.2 =0,    (9)

    ∂.sup.3 x.sub.r /∂t.sup.3 =ω.sub.e.sup.3 {x.sub.r0 sin ω.sub.e (t-t0)-y.sub.r0 cos ω.sub.e (t-t0)},(10)

    ∂.sup.3 y.sub.r /∂t.sup.3 =ω.sub.e.sup.3 {x.sub.r0 cos ω.sub.e (t-t0)+y.sub.r0 sin ω.sub.e (t-t0)},(11)

    ∂.sup.3 z.sub.r /∂t.sup.3 =0.    (12)

The LD signal source(s) is/are assumed to be moving, with known location coordinates as a function of time. For definiteness, each LD signal source is assumed to be carried on a satellite that revolves around the Earth in a circular orbit of radius a with an orbital plane that includes the center of the Earth. Each satellite rotates with a known characteristic angular velocity ω_(s), in an orbital plane that is described by a known azimuthal angle φ and by a known angle θ that this plane makes with the equatorial line on the Earth. A satellite thus has a time-varying location given by the Cartesian coordinates

    x.sub.s (t)=a cos φ cos θ cos ω.sub.s t"+a sin φ sin ω.sub.s t",                                         (13)

    y.sub.s (t)=-a sin φ cos θ cos ω.sub.s t"+a cos φ sin ω.sub.s t",                                         (14)

    z.sub.s (t)=a sin θ cos ω.sub.s t".            (15)

The time variable t" in Eqs. (13)-(15) is a translated time variable t"=t-t(j), to account for the difference in location of each satellite (j). The parameters a, ω_(s), φ and θ also vary with the satellite (j). The first, second and third time derivatives of the location coordinates of the satellite are found to be

    ∂x.sub.s /∂t=aω.sub.s {-cos φ cos θ sin ω.sub.s t"+sin φ cos ω.sub.s t"},(16)

    ∂y.sub.s /∂t=aω.sub.s {-sin φ cos θ sin ω.sub.s t"+cos φ cos ω.sub.s t"},(17)

    ∂z.sub.s /∂t=-aω.sub.s sin θ sin ω.sub.s t",                                         (18)

    ∂.sup.2 x.sub.s /∂t.sup.2 =aω.sub.s.sup.2 {-cos φ cos θ cos ω.sub.s t"-sin φ sin ω.sub.s t"},                                                      (19)

    ∂.sup.2 y.sub.s /∂t.sup.2 =aω.sub.s.sup.2 {-sin φ cos θ cos ω.sub.s t"-cos φ sin ω.sub.s t"},                                                      (20)

    ∂.sup.2 z.sub.s /∂t.sup.2 =-aω.sub.s.sup.2 sin θ cos ω.sub.s t",                         (21)

    ∂.sup.3 x.sub.s /∂t.sup.3 =aω.sub.s.sup.3 {cos φ cos θ sin ω.sub.s t"-sin φ cos ω.sub.s t"},                                                      (22)

    ∂.sup.3 y.sub.s /∂t.sup.3 =aω.sub.s.sup.3 {-sin φ cos θ sin ω.sub.s t"-cos φ cos ω.sub.s t"},                                                      (23)

    ∂.sup.3 z.sub.s /∂t.sup.3 =-aω.sub.s.sup.3 sin θ sin ω.sub.s t".                         (24)

A computer 14 is programmed to handle the computations.

The square of the measured pseudorange is first formed. If the vectors r_(r) (t) and r_(s) (t) represent the respective vectors from the coordinate origin to the LD signal receiver (i) and to the satellite (j), the (corrected) pseudorange equation for a selected receiver (i) and a selected satellite (j) becomes ##EQU1## The unknowns in Eq. (25) are the location coordinates x_(r0), y_(r0), z_(r0) and the receiver time offset b_(r) (optionally determined). After simplification, the time derivative of the separation distance (squared) becomes ##EQU2##

    m.sub.x (t)=a{ω.sub.e cos θ sin (ω.sub.e t'-φ) cos ω.sub.s t"+ω.sub.e cos (ω.sub.e t'-φ)

    sin ω.sub.s t"+ω.sub.s cos θ cos (ω.sub.e t'-φ) sin ω.sub.s t"+ω.sub.s sin (ω.sub.e t'-φ) cos ω.sub.s t"},                                        (27)

    m.sub.y (t)=a{-ω.sub.e cos θ cos (ω.sub.e t'-φ) cos ω.sub.s t"+ω.sub.e sin (ω.sub.e t'-φ)

    sin ω.sub.s t"+ω.sub.s cos θ sin (ω.sub.e t'-φ) sin ω.sub.s t"-ω.sub.s cos (ω.sub.e t'-φ) cos ω.sub.s t"},                                        (28)

    m.sub.z (t)=ω.sub.s sin θ sin ω.sub.s t, (29)

    Δv=∂r.sub.r /∂t-∂r.sub.s /∂t,                                         (30)

More generally, an nth power of the quantity (Δr·Δv) becomes

    (Δr·Δv).sup.n ={m.sub.x (t) x.sub.r0 +m.sub.y (t) y.sub.r0 +m.sub.z (t) z.sub.r0 }.sup.n.                   (31)

The first time derivative of Doppler shift involves the quantity ##EQU3##

    Δu=∂.sup.2 r.sub.r /∂t.sup.2 -∂.sup.2 r.sub.s /∂t.sup.2.     (33)

The second derivative of Doppler shift involves the quantity ##EQU4##

    Δw=∂.sup.3 r.sub.r /∂t.sup.3 -∂.sup.3 r.sub.s /∂t.sup.3.     (35)

For a first body whose motion can be described by the usual special relativistic equations, the frequency Doppler shift of a signal emitted by that body and received by a second body is given by

    Δν(Doppler;t;i;j)=-ν.sub.0 (Δr·Δv)/c|Δr|,(36)

(C. Moller, The Theory of Relativity, Oxford Univ. Press, 1952, pp. 58-62), where Δr and Δv are the first body-second body location vector difference and velocity vector difference, respectively, and ν₀ is the signal frequency as transmitted by the first body.

Using the cosine law, the quantity

    |Δr|=|r.sub.r (t;i)-r.sub.s (t;j)|(37)

in the denominator in Eq. (45) can be replaced by an approximation

    |Δr|={a.sup.2 +a.sub.e.sup.2 -2 a a.sub.e cos ψ}.sup.1/2,                                           (38)

where a_(e) is the Earth's radius and ψ is the estimated angle between a vector directed toward the LD signal receiver and a vector directed toward the satellite whose signal is being received (π/2≦ψ≦π). The numerical value |Δr| may be written as χ(a-a_(e)), where a=26,560 km, a_(e) =7,040 km and 1≦χ≦χ(max)={(a+a_(e))/(a-a_(e))}^(1/2) =1.312 for GPS orbits, and a=25,510 km, a_(e) =7,040 km and 1≦χ≦χ(max)=1.328 for GLONASS orbits, for example. Preferably, the quantity |Δr| is estimated using a measured pseudorange value.

The first and second time derivatives of the frequency Doppler shift between the first body and the second body are given by

    ∂{c|Δr|Δν/ν0}/.differential.t=(c|Δr|/ν.sub.0) ∂Δν/∂t-(cΔν/ν.sub.0).sup.2 =-Δv·Δv-Δr·Δu,  (39)

    ∂.sup.2 {c|Δr|Δν/ν0}/∂t.sup.2 =(c|Δr|/ν.sub.0)∂.sup.2 Δν/∂t.sup.2 -3(c/ν.sub.0).sup.2 Δν∂Δν/∂t=-2Δv·.DELTA.u-Δr·Δw.                        (40)

In a first embodiment of the invention, a pseudorange (optional) and Doppler shifts Δν(Doppler;t;i;j) are measured for three different (non-collinear) satellites, j=1, 2, 3, at a selected time t=t1. From Eqs. (26), the three Doppler shift relations become ##EQU5## where the quantities in the 3×3 matrix M and in the 3×1 column vector G are calculated or measured. Equation (40) is then inverted to provide the receiver (initial) location coordinates ##EQU6## The matrix M(t1;i;1,2,3) is invertible, in part because the three satellites are non-collinear. Equations (1), (2) and (3) can be applied to estimate the receiver location coordinates on the rotating Earth at any other time t. If the receiver clock offset b_(r) is also needed (optional), the pseudorange PR(t1;i;j) is also measured, and Eq. (25) is applied to determine the offset value br.

In a second embodiment, three or four of pseudorange and the zeroth, first and second time derivatives of the Doppler shift for a single satellite (j) are measured and/or calculated for a selected time t=t1. The three Eqs. (36), (39) and (40), supplemented by Eqs. (26), (32) and (34), can be re-expressed in the forms ##EQU7## where M'(t1;i;j)=(M'_(k1)) is a 3×3 matrix with entries

    M'.sub.11 =m.sub.x (t1;i;j),                               (48)

    M'.sub.12 =m.sub.y (t1;i;j),                               (49)

    M'.sub.13 =m.sub.z (t1;i;j),                               (50)

    M'.sub.21 =∂m.sub.x (t1;i;j)/∂t, (51)

    M'.sub.22 =∂m.sub.y (t1;i;j)/∂t, (52)

    M'.sub.23 =∂m.sub.z (t1;i;j)/t,               (53)

    M'.sub.31 =∂.sup.2 m.sub.x (t1;i;j)/∂t.sup.2,(54)

    M'.sub.32 =∂.sup.2 m.sub.y (t1;i;j)/∂t.sup.2,(55)

    M'.sub.33 =∂.sup.2 m.sub.z (t1;i;j)/∂t.sup.2,(56)

Equation (46) is then inverted to provide the receiver (initial) location coordinates ##EQU8## The matrix M'(t1;i;1,2,3) should be verified to be invertible by computation of its determinant. Equations (1), (2) and (3) can be applied to estimate the present receiver location coordinates on the rotating Earth at any time t. Equation (25) is applied to estimate the time offset b_(r) for the receiver, if this quantity is needed. Equations (41), (46) and (71)-(80) are exact if separation distance |Δr(t1;i;j)| can be measured or otherwise determined. If receiver clock offset b_(r) is also needed, the pseudorange PR(t1;i;j) is also measured or calculated, and Eq. (25) is applied to determine the offset value b_(r).

In a third embodiment of the invention, pseudorange PR(t;i;j) is measured, preferably for at least two satellites, j=j1 and j=j2, and Eq. (25) is rewritten for each of these two satellites in the forms

    (x.sub.r (t;i)-x.sub.s (t;j1)).sup.2 +(y.sub.r (t;i)-y.sub.s (t;j1)).sup.2 +(z.sub.r (t;i)-z.sub.s (t;j1)).sup.2 =(PR(t;i;j1)-b.sub.r).sup.2,(58)

    (x.sub.r (t;i)-x.sub.s (t;j2)).sup.2 +(y.sub.r (t;i)-y.sub.s (t;j2)).sup.2 +(z.sub.r (t;i)-z.sub.s (t;j2)).sup.2 =(PR(t;i;j2)-b.sub.r).sup.2.(59)

Equation (58) is subtracted from Eq. (59), and the result is rearranged to yield ##EQU9##

    Δx.sub.s (t;j1;j2)=x.sub.s (t,j1)-x.sub.s (t;j2),    (61)

    Δy.sub.s (t;j1;j2)=y.sub.s (t;j1)-y.sub.s (t;j2),    (62)

    Δz.sub.s (t;j1;j2)=z.sub.s (t;j1)-z.sub.s (t,j2),    (63)

    ΔPR(t;i;j1;j2)=PR(t;i;j1)-PR(t;i;j2),                (64)

    r.sub.s (t;j)2=x.sub.s (t;j).sup.2 +y.sub.s (t;j).sup.2 +z.sub.s (t;j).sup.2.                                              (65)

Equation (60) is rewritten in terms of the unknown quantities x_(r0), y_(r0) and z_(r0) as ##EQU10##

    Q(t;i;j1;j2)={r.sub.s (t,j1).sup.2 -r.sub.s (t;j2).sup.2 +PR(t;i;j2).sup.2 -PR(t;i;j1).sup.2 }/2,                                    (67)

At least two of the following six quantities are also measured: the zeroth, first and second time derivatives of Doppler shift for satellite j=j1;

and the zeroth, first and second time derivatives of the Doppler shift for j=j2.

From Eqs. (26) and (46), the zeroth derivative of Doppler shift for a signal received from satellite j may be rearranged to

    m.sub.x (t;i;j)x.sub.r0 +m.sub.y (t;i;j)y.sub.r0 +m.sub.z (t;i;j)z.sub.r0 =-Δνc|Δr(t;i;j)|/ν.sub.0,(68)

where the quantities m_(x) (t), m_(y) (t), m_(z) (t) are set forth in Eqs. (27)-(29). From Eq. (39), a first selected combination of zeroth and first time derivatives of the Doppler shift may be expressed as

    (∂m.sub.x (t;i;j)/∂t)x.sub.r0 +(∂m.sub.y (t;i;j)/∂t)y.sub.r0 +(∂m.sub.z (t;i;j)/∂t)z.sub.r0 =-(c|Δr(t1;i;j)|/ν.sub.0)(∂Δν/∂t)+(cΔν/ν.sub.0).sup.2.     (69)

From Eqs. (40), a second selected combination of zeroth, first and second time derivatives of the Doppler shift may be expressed as

    (∂.sup.2 m.sub.x (t;i;j)/∂t.sup.2)x.sub.r0 +(∂.sup.2 m.sub.y (t;i;j)/∂t.sup.2)y.sub.r0 +(∂.sup.2 m.sub.z (t;i;j)/∂t.sup.2)z.sub.r0=-(c|Δr(t1;i;j)|/ν.sub.0)(∂.sup.2 Δν/∂t.sup.2)+(3c.sup.2 /ν.sub.0.sup.2)Δν(∂Δν/∂t).(70)

Each of Eqs. (68), (69) and (70) can be expressed in the general form

    M.sub.x,k (t;i;j)x.sub.r0 +M.sub.y,k (t;i;j)y.sub.r0 +M.sub.z,k (t;i;j)z.sub.r0 =G.sub.k (t;i;j),                         (71)

where the index k=0, 1 and 2 corresponds, respectively, to the zeroth time derivative of Doppler shift, the first selected combination of zeroth and first time derivatives of Doppler shift, and the second selected combination of zeroth, first and second time derivatives of Doppler shift, respectively, for a signal received at the receiver (i) from the satellite (j). Because Doppler shift measurements from two satellites (j=j1 and j=j2) are available here, Eq. (71) represents any of six Doppler shift measurement combinations: three for the satellite j=j1 and three for the satellite j=j2.

Two of the six relations (71) are chosen for index values k=k1 and k=k2 and for the satellite indices j=j1 and j=j2, respectively, where j1 may be the same as, or may be different from, j2. If j1=j2, the indices k1 and k2 must differ; if j1≠j2, the index k1 may be the same as, or may be different from, the index k2. These two relations are rewritten in the form ##EQU11## and Eq. (72) is inverted to express x_(r0) and y_(r0) in terms of the unknown z_(r0) :

    x.sub.r0 ={M.sub.y,k2 (t;i;j2)G.sub.k1 (t;i;j1)--M.sub.y,k1 (t;i;j1)G.sub.k2 (t;i;j2)-M.sub.z,k1 (t;i;j2)z.sub.r0 }/(det1),(73)

    y.sub.r0 ={-M.sub.x,k2 (t;i;j2)G.sub.k1 (t;i;j1)-M.sub.z,k1 (t;i;j1)z.sub.r0 +M.sub.x,k1 (t;i;j1)G.sub.k2 (t;i;j2)-M.sub.z,k1 (t;i;j2)z.sub.r0 }/(det1),                                (74)

    det1=M.sub.x,k1 (t;i;j1)M.sub.y,k2 (t;i;j2)-M.sub.x,k2 (t;i;j2)M.sub.y,k1 (t;i;j1).                                                 (75)

For simplicity in notation, Eqs. (73) and (74) are rewritten as

    x.sub.r0 =C.sub.x (t;i;j1;j2)+D.sub.x (t;i;j1;j2)z.sub.r0, (76)

    y.sub.r0 =C.sub.y (t;i;j1;j2)+D.sub.y (t;i;j1;j2)z.sub.r0, (77)

with the coefficients C_(x), D_(x), C_(y) and D_(y) thus defined.

Equations (76) and (77) are now used in Eq. (66) to express the unknown z_(r0) in terms of the remaining unknown, b_(r), ##EQU12## For simplicity of notation, Eq. (77) is rewritten as

    z.sub.r0 =E.sub.z (t;i;j1;j2)+F.sub.z (t;i;j1;j2)b.sub.r.  (79)

Equations (76), (77) and (79) are then substituted in Eq. (25) to produce a quadratic equation ##EQU13## for the remaining unknown, the clock offset b_(r), where t'=t-t0 and t"=t-t(j). The quadratic relation (80) is then solved for the clock offset b_(r) ; only the more realistic value of b_(r), close to zero, is accepted here. This value of b_(r) is used in Eq. (79) to determine the value of z_(r0), and the value for z_(r0) is used in Eqs. (76) and (77) to obtain the values for x_(r0) and y_(r0). Note that the value for clock offset b_(r) is, or can be, determined in this approach for the third embodiment, whether or not this value is of any intrinsic interest. In other approaches for this third embodiment, a quadratic equation can be developed in terms of one of the other unknowns, x_(r0), y_(r0) or z_(r0), using a similar development, and the value of the clock offset b_(r) may not need to be computed explicitly.

To recapitulate, this third embodiment uses two pseudorange measurements, for satellites j=j1 and j=j2, and uses two of the following six combinations of measurements: (1) the zeroth derivative of Doppler shift for satellite j=j1 and/or for satellite j=j2; (2) a first selected combination of the zeroth and first time derivatives of Doppler shift for satellite j=j1 and/or for satellite j=j2; and (3) a second selected combination of the zeroth, first and second time derivatives of Doppler shift for satellite j=j1 and/or for satellite j=j2, where the index values j1 and j2 are chosen independently. An alternate version uses one pseudorange value to obtain b_(r) (optional) and three of the six combinations to obtain x_(r0), y_(r0) and z_(r0).

Thus far, the development has focused primarily on measurement and analysis of Doppler shift and its time derivatives. This development is extended to pseudorange and its time derivatives by the following observations. Pseudorange PR, after removal of error terms such as satellite clock offset, satellite noise, receiver noise, ionospheric propagation time delay, tropospheric propagation time delay and multipath errors, is approximately equal to the magnitude of the radius vector difference,

    PR(t;i;j)≈=|Δr(t;i;j)|=|r.sub.r (t;i)-r.sub.s (t;j)|.                            (81)

The first time derivative of pseudorange becomes

    ∂|Δr|/∂t=Δr.multidot.Δv/|Δr|=-(c/ν.sub.0)Δν,(82)

which is the Doppler shift, apart from a constant of proportionality, -c/ν₀.

Thus, measurement of the zeroth, first and second time derivatives of Doppler shift is substantially the same as measurement of the first, second and third time derivatives, respectively, of pseudorange. One can thus work with multiple time derivatives of pseudorange and/or with multiple time derivatives of Doppler shift (including the zeroth derivative) within the same formalism. For ease of reference, pseudorange measurement values and Doppler shift measurement values from a selected satellite to an LD signal receiver are collectively referred to as "location signal measurement" (LSM) values.

FIG. 2 is a flow chart illustrating the first embodiment of the invention. In step 21, at least one LSM value (pseudorange or Doppler shift) is determined for three non-collinear satellites. In step 23 (optional), LSM values for the three satellites (e.g., entries of the matrix M(t1;i;1;2;3) in Eq. (42)) are determined. In step 25, the initial location coordinates x_(r0), y_(r0) and z_(r0) are determined, using Eq. (45) or equivalent. In step 27 (optional), the receiver time offset b_(r) is determined, using the pseudorange measurement in Eq. (25) or some other approach. In step 29 (optional), the present location coordinates x_(r) (t), y_(r) (t) and z_(r) (t) are determined, using Eqs. (1), (2) and (3) or some other suitable equations.

FIG. 3 is a flow chart illustrating the second embodiment. In step 31, three or four LSM and/or LSM derivative values (pseudorange, zeroth, first and second time derivatives of Doppler shift) are determined for a single satellite. In step 33, LSM and derivative values for the entries of the matrix M'(t1;i;j) in Eq. (46) are used to determine the initial location coordinates x_(r0), y_(r0) and z_(r0), using Eqs. (46) or (57). In step 35 (optional), the receiver time offset b_(r) is determined, using Eq. (25) or some other approach. In step 37 (optional), the present location coordinates x_(r) (t), y_(r) (t) and z_(r) (t) are determined, using Eqs. (1), (2) and (3) or some other suitable equations.

FIG. 4 is a flow chart illustrating the third embodiment of the invention. In step 41, a first LSM value (preferably pseudorange) for each of two satellites (j=j1 and j=j2) is measured. In step 43, a first linear relation (corresponding to Eq. (66)) is obtained between the four unknowns x_(r0), y_(r0), z_(r0), and b_(r). In step 45, third and fourth LSM values (two of the six Doppler shift time derivative combinations) are measured, and second and third linear relations are obtained between three of the unknowns x_(r0), y_(r0) and z_(r0). In step 47, a quadratic relation involving only one of the four unknowns x_(r0), y_(r0), z_(r0), and b_(r) is obtained, using the first, second and third linear relations, and a (realistic) solution is determined. In step 49, the remaining initial receiver location fix coordinates, x_(r0), y_(r0), z_(r0), and b_(r), (or of three location fix coordinates, if the value of b_(r) or of one of the other location fix coordinates is not of interest) are computed. In step 51 (optional), the present location coordinates x_(r) (t), y_(r) (t) and z_(r) (t) are determined, using Eqs. (1), (2) and (3) or some other suitable equations

For any embodiment, if three or fewer receiver location fix coordinates are to be estimated, for example, because the remaining coordinates are known, the total number of zeroth, first and second time derivatives of Doppler shift is reduced correspondingly. For example, if only two location coordinates are estimated, only two of the three Doppler shift derivatives need be measured from a single satellite, or from one Doppler shift derivative from each of two satellites. Measurement of pseudorange from one satellite plus Doppler shift(s) allows determination of clock offset. If pseudoranges are measured from two or more satellites, the number of Doppler shift derivative measurements is reduced by two, as shown in the preceding development.

The time derivatives for an LSM value can be estimated by at least two approaches. In a first, straightforward approach, LSM values L(t;i,j) are measured from signals received at a receiver (i) from a selected LD signal source (j) at a sequence of times t=τ1, τ2, τ3, τ4, . . . that satisfy τ1<τ2≦τ3<τ4< . . . A first time derivative of the LSM value at t=τ2 is defined as ##EQU14## where the dependence upon receiver index (i) and upon LD signal source index (j) are suppressed here to simplify the notation. A second time derivative of LSM value at a time t=τ2'=(τ2+t3)/2 is defined as ##EQU15## Users of this first approach to determining LSM value time derivatives may experience difficulties with reduced signal-to-noise (S/N) ratios because of the definitions adopted for the derivatives in Eqs. (83) and (84).

In a second approach, the S/N ratio is not reduced as severely, but the required computation is increased. In a first variation of the second approach, the LSM value is measured at several distinct values t=τk(k=1, 2, 3, 4,. . . ), with τ1<τ2<τ3<τ4<. . . , and LSM values for intermediate values (τ1<t<τ4) are defined by interpolation as

    L(t;2)=L(τ2)+A(t-τ2)+B(t-τ2).sup.2 +C(t-t2).sup.3 (τ1≦t≦τ4),                          (85) ##EQU16##

    det3=(τ3-τ1)(τ3-τ2)(τ2-τ1)(τ4-τ2)(τ4-.tau.3)(τ3-τ2),                                      (87)

and the time derivatives at t=τ2 become

    dL(t=τ2)/dt=A,                                         (88)

    d.sup.2 L(t=τ2)/dt.sup.2 =2B,                          (89)

    d.sup.3 L(t=τ2)/dt.sup.3 =6C.                          (90)

Polynomials of degree higher than the second, third or fourth degree polynomials used in Eqs. (85) and (92) can also be used to determine the first and second derivatives of the LSM values, by an obvious extension of the development in Eqs. (85)-(90).

In a third approach, M not-necessarily-consecutive values (M≧3) of LSM values L(τm) are measured and a polynomial in t of degree M'(≦M-2) is chosen to pass as closely as possible, according to a selected criterion, to each of the measured LSM values. For example, the polynomial may be ##EQU17## and the criterion may be least squares minimization of the collective error associated with the M measured values, by suitable choices of the polynomial coefficients p_(k) (k=0, 1, . . . , M'): ##EQU18## Equations (93) provide a linear system of M'+1 relations in the polynomial coefficients p_(j), which can be solved by standard linear algebra techniques. The first and second derivatives of LSM values are then estimated as ##EQU19## which can be evaluated for any selected time t (with arbitrary time origin).

Whether the first approach, the second approach, the third approach, or some other reasonable approach is used here to estimate the first and second derivatives of LSM values from (only) the measured LSM values, the estimates for these derivatives will each be a linear combination of a plurality of measured values L(τk) (k=1, 2, 3, 4, . . . ). At least two such measured values and at least three measured values of LSM values are required to estimate the first derivative and the second derivative, respectively.

Another method for obtaining equations relating LSM values and their time derivatives to (stationary) receiver location coordinates makes use of the circular rotation of the receiver about the Earth's rotational axis and of the circular rotation of the satellite (LD signal source) about the Earth. If r_(r) (t;i) and r_(s) (t;j) represent the present radius vectors (relative to the Earth's center or another common origin) of the receiver (i) and the satellite (j), respectively, the time derivatives of these quantities may be expressed as

    dr.sub.r (t;i)/dt=ω.sub.r r.sub.r (t;i),             (97)

    d.sup.n r.sub.r (t;i)/dt.sup.n =ω.sub.r d.sup.n-1 r.sub.r (t;i)/dt.sup.n-1 (n=2, 3, . . .),                         (98)

    dr.sub.s (t;j)/dt=ω.sub.s r.sub.s (t;j),             (99)

    d.sup.n r.sub.s (t;j)/dt.sup.n =ω.sub.s d.sup.n-1 r.sub.s (t;j)/dt.sup.n-1 (n=2, 3, . . .),                         (100)

where ω_(r) or and ω_(s) are the respective angular velocity vectors (constant, with dimensions sec-⁻¹) associated with the rotation of the (stationary) receiver about the Earth's axis and with the rotation of the satellite about the Earth. Using the notations introduced in Eqs. (25), (30), (33) and (35), the zeroth, first and second Doppler shift derivatives are expressed as ##EQU20##

Each of the Eqs. (101)-(103) is a scalar equation, linear in the three unknown components (x_(r) (t;i),y_(r) (t;i),z_(r) (t;i)) of the vector r_(r), with coefficients that are either constant or are known functions of t. The left hand sides of the Eqs. (101)-(103) are measurable combinations of the zeroth, first and second derivatives of Doppler shift (multiplied in some instances by the quantity |Δr(t;i;j)|). These three scalar equations are expressed in matrix form as

    M"(t;i;j1;j2;j3)X03=G"(t;i,j1;j2;j3),                      (104)

    M".sub.11 =(ω.sub.s,y -ω.sub.r,y)z.sub.s (t,j1)-(ω.sub.s,z -ω.sub.r,z)y.sub.s (t;j1),  (105)

    M".sub.12 =(ω.sub.s,z -ω.sub.r,z)x.sub.s (t,j1)-(ω.sub.s,x-ω.sub.r,x)z.sub.s (t;j1),   (106)

    M".sub.13 =(ω.sub.s,x -ω.sub.r,x)y.sub.s (t;j1)-(ω.sub.s,y-ω.sub.r,y)x.sub.s (t;j1),   (107)

    M".sub.21 =(ω.sub.r.sup.2 +ω.sub.s.sup.2)x.sub.s (t;j2)+(ω.sub.r ·r.sub.s) ω.sub.r,x +(ω.sub.s ·r.sub.s)ω.sub.s,x,                        (108)

    M".sub.22 =(ω.sub.r.sup.2 +ω.sub.s.sup.2)y.sub.s (t;j2)+(ω.sub.r ·r.sub.s) ω.sub.r,y +(ω.sub.s ·r.sub.s)ω.sub.s,y,                        (109)

    M".sub.23 =(ω.sub.r.sup.2 +ω.sub.s.sup.2)z.sub.s (t;j2)+(ω.sub.r ·r.sub.s) ω.sub.r,z +(ω.sub.s ·r.sub.s)ω.sub.s,z,                        (110) ##EQU21##

    G".sup.tr ={G".sub.1, G".sub.2, G".sub.3 }.sup.tr,         (114)

    G".sub.1 =c|Δr|Δν(t;i;j1)/ν.sub.0,(115)

    G".sub.2 =(c/ν.sub.0)|Δr |dΔν(t;i;j2)/dt-(cΔν(t;i;j2)/ν.sub.0).sup.2,(116)

    G".sub.3 =(c|Δr|/ν.sub.0)d.sup.2 Δν(t;i,j3)/dt.sup.2 -3(c/ν.sub.0).sup.2 Δν(t;i,j3)dΔν/dt,                       (117) ##EQU22## Here, the satellite indices j1, j2 and j3 may be the same or may differ from each other, indicating that Doppler shift derivative measurements may be made using timed signals from the same satellite or from different satellites. The solution matrix X3 is found, as usual, by inverting the matrix M(t;i;j1; j2; j3):

    X3=M(t;i;j1;j2;j3).sup.-1 G"(t;i;j1;j2;j3).                (119)

Here, the solutions are the present values (x_(r) (t;i),y_(r) (t;i),z_(r) (t;i)) of the three location coordinates, not the "initial" values (x_(r) (t0;i),y_(r) (t0;i),z_(r) (t0;i)) of these coordinates. Note that no restrictions are imposed on the angular velocity vectors ω_(r) and ω_(s), except that each is referenced to the same origin, such as the Earth's center.

The orbit(s) of the satellite(s), assumed to be circular in Eqs. (13)-(15), may be extended to an elliptical orbit in the first embodiment, illustrated in FIG. 5, by the following development. The ellipse has a semi-major axis (x') of length a, a semi-minor axis (parallel to y') of length b, and an eccentricity ε given by

    ε={1-b.sup.2 /a.sup.2 }.sup.1/2,                   (120)

where b<a. The radius length ρ of a line segment from an ellipse focus F0 to its intersection with the ellipse curve is given by

    ρ(φ)=a(1-ε.sup.2)/(1+ε cos (φ+π/2)),(121)

where φ is the angle the line segment makes with the y'-axis in FIG. 5. (K. R. Symon, Mechanics, Addison Wesley, Cambridge, Mass., 1953, pp. 111-114). Note that the azimuthal angle φ is defined with reference to the y'-axis rather than to the more conventional x'-axis. The analogues of Eqs. (13)-(15) for an elliptical orbit, instead of a circular orbit, are verified to be

    x.sub.s (t;j)=ρ(φ) cos Φ cos θ cos ω.sub.s (t-t(j))+ρ(φ) sin Φ sin ω.sub.s (t-t(j)),(122)

    y.sub.s (t;j)=ρ(φ) sin Φ cos θ cos ω.sub.s (t-t(j))+ρ(φ) cos Φ sin ω.sub.s (t-t(j)),(123)

    z.sub.s (t;j)=ρ(φ) sin θ cos ω(t-t(j)),(124)

    ρ(φ).sup.2 (dφ/dt)=constant (κ),         (125)

where the "anchoring" body is located at an ellipse focus FO (FIG. 5) and κ is a measure of the angular momentum of the orbiting satellite. Using a result found in G. Petit Bois, Tables of Indefinite Integrals, Dover Publs. New York, 1961, p. 118, Eq. (131) can be integrated to yield an implicit relation between φ and t for an elliptic orbit:

    κt=-ε cos φ/{(1-ε.sup.2)(1-ε sin φ)}+2(1-ε.sup.2).sup.-3/2  tan.sup.-1 {(tan(φ/2)-ε)/(1-ε.sup.2).sup.1/2 }.  (126)

In some situations, it may not be possible to obtain sufficiently accurate measurements on signals from one or more satellites to provide values for pseudorange, Doppler shift and two time derivatives of Doppler shift (or, equivalently, of pseudorange and three time derivatives of pseudorange). For example, it may be possible only to obtain values of pseudorange and the first and second time derivatives of pseudorange for a satellite at the receiver. In this situation, receiver location can be estimated, using an approximation. Alternatively, this approximation may be used if one of the location fix coordinates, receiver clock offset b_(r), is known to be small and is not of much interest. In this approximation, one returns to Eq. (25) and re-expresses this as an approximation, ##EQU23## where presence of the receiver clock offset b_(r) is ignored. The linear relations (73) are then used to express the quantities x_(r0) and y_(r0) in terms of z_(r0), and Eq. (133) is restated as a quadratic equation in the remaining unknown, here z_(r0). ##EQU24## The more realistic solution z_(r0) for Eq. (134) is then accepted and used in Eqs. (73)-(74) to determine the remaining initial location coordinates x_(r0) and y_(r0). Equations (1)-(3) can be used, as before, to determine the values of the x-, y- and z-coordinates of the receiver at some other selected time t.

FIG. 6 is a schematic view of apparatus suitable for practicing the invention and may be a console or be handheld. The apparatus includes an LD signal antenna 12 that receives and feeds LD signals to an LD signal receiver 11 for processing and determination of the LSM values. The LSM values are then passed to a computer 14, which may be part of the LD signal receiver 11 or may operate independently. The computer 14 includes a microprocessor 15 and associated memory module 16 that are programmed to perform the computations on the pseudorange values and Doppler shift time derivative values to determine the initial receiver location fix coordinate values (x_(r0),y_(r0) z_(r0),b_(r)), receiver coordinate values (x_(r) (t),y_(r) (t),z_(r) (t)) at some other selected time, elliptical orbit parameters, or any other relevant numerical quantity. The computer 14 optionally includes a data entry module 17, such as a keyboard or light pen, a visual and/or audible display 18 for the information computed or measured, and a power module to provide power for the other components shown in FIG. 6. 

We claim:
 1. A method for estimating location coordinates of a signal receiver from signals received from three satellites, the method comprising the steps of:receiving a timed signal from each of first, second and third satellites at a receiver at approximately a first selected time; determining at least one of a pseudorange value and a Doppler shift value, with these values being collectively denoted location signal measurement values or LSM values, for the timed signal received from each of the three satellites, where at least one of the LSM values is a Doppler shift value; and using the determined LSM values to estimate first and second location coordinate values for the receiver and at least one of a receiver time offset value and a third receiver location coordinate value, at a second selected time, determined with reference to the first selected time.
 2. The method of claim 1, further comprising the steps of:determining a selected LSM value for said timed signal for a selected one of said satellites at approximately said first selected time, where at least one of said first, second and third LSM values and the selected LSM value is a pseudorange value; and estimating said receiver time offset value and said third location coordinate value from the selected LSM value and said estimated first and second location coordinate values for said receiver.
 3. A method for estimating location coordinates of a signal receiver from signals received from one satellite, the method comprising the steps of:receiving a timed signal from a selected satellite at a receiver at approximately a first selected time; determining at least two of a pseudorange value, a Doppler shift value, a first selected combination of Doppler shift value and first time derivative of Doppler shift value, and a second selected combination of Doppler shift value and first and second time derivatives of Doppler shift value, with these values being collectively denoted location signal measurement values or LSM values, for the timed signal received from the satellite; using the at least two determined LSM values to estimate first and second location coordinate values for the receiver, at a second selected time determined with reference to the first selected time.
 4. The method of claim 3, further comprising the steps of:determining a third LSM value, for the timed signal received from the satellite, at approximately said first selected time; and using the determined third LSM value and said estimated first and second location coordinate values to estimate at least one of a receiver time offset value and a third location coordinate value for said receiver, at said second selected time.
 5. The method of claim 3, further comprising the steps of:determining third and fourth LSM values, for the timed signal received from the satellite, at approximately said first selected time, where at least one of said at least two LSM values and the third and fourth LSM values is a pseudorange value; and using the third and fourth LSM values and said estimated first and second location coordinate values to estimate a receiver time offset value and a third location coordinate value for said receiver, at said second selected time.
 6. A method for estimating location coordinates of a signal receiver from signals received from two satellites, the method comprising the steps of:receiving a timed signal from each of first and second satellites at a receiver at approximately a first selected time; determining at least one of pseudorange value, Doppler shift value, a first selected combination of Doppler shift value and first time derivative of Doppler shift value, and a second selected combination of Doppler shift value and first and second time derivatives of Doppler shift value, with these values being collectively denoted location signal measurement values or LSM values, for the timed signal received from the first satellite; determining an LSM value for the timed signal received from the second satellite; and using the LSM values from the satellites to estimate first and second location coordinate values for the receiver, at a second selected time, determined with reference to the first selected time.
 7. The method of claim 6, further comprising the steps of:determining a selected LSM value for said timed signal for a selected one of said satellites at approximately said first selected time; and estimating at least one of a receiver time offset value and a third receiver location coordinate value from the selected LSM value and said estimated first and second location coordinate values for said receiver.
 8. The method of claim 6, further comprising the steps of:determining a second LSM value for said timed signal for each of said first and second satellites at approximately said first selected time, where at least one of said first LSM values and the second LSM values for said two satellites is a pseudorange value; and estimating a receiver time offset value and a third location coordinate value for said receiver from the second LSM values and said estimated first and second location coordinate values for the receiver.
 9. Apparatus for estimating location coordinates of a signal receiver from signals received from three satellites, the apparatus comprising:a location determination unit that: receives a timed signal from each of first, second and third satellites at a receiver at approximately a first selected time; and determines at least one of a pseudorange value and a Doppler shift value, with these values being collectively denoted location signal measurement values or LSM values, for the timed signal received from each of the three satellites, where at least one of the LSM values is a Doppler shift value; and a computer that is programmed to use the determined LSM values to estimate first and second location coordinate values for the receiver and at least one of a receiver time offset value and a third receiver location coordinate value, at a second selected time, determined with reference to the first selected time.
 10. The apparatus of claim 9, wherein:said location determination unit determines a selected LSM value for said timed signal for a selected one of said satellites at approximately said first selected time, where at least one of said first, second and third LSM values and the selected LSM value is a pseudorange value; and said computer is further programmed to estimate said receiver time offset value and said third location coordinate value from the selected LSM value and said estimated first and second location coordinate values for said receiver.
 11. Apparatus for estimating location coordinates of a signal receiver from signals received from one satellite, the apparatus comprising:a location determination unit that:receives a timed signal from a selected satellite at a receiver at approximately a first selected time; and determines at least two of a pseudorange value, a Doppler shift value, a first selected combination of Doppler shift value and first time derivative of Doppler shift value, and a second selected combination of Doppler shift value and first and second time derivatives of Doppler shift value, with these values being collectively denoted location signal measurement values or LSM values, for the timed signal received from the satellite; a computer that is programmed to use the at least two determined LSM values to estimate first and second location coordinate values for the receiver, at a second selected time determined with reference to the first selected time.
 12. The apparatus of claim 11, wherein:said location determination unit determines a third LSM value, for the timed signal received from the satellite, at approximately said first selected time; and said computer is further programmed to use the determined third LSM value and said estimated first and second location coordinate values to estimate at least one of a receiver time offset value and a third location coordinate value for said receiver, at said second selected time.
 13. The apparatus of claim 11, wherein:said location determination unit determines third and fourth LSM values, for the timed signal received from the satellite, at approximately said first selected time, where at least one of said at least two LSM values and the third and fourth LSM values is a pseudorange value; and said computer is further programmed to use the third and fourth LSM values and said estimated first and second location coordinate values to estimate a receiver time offset value and a third location coordinate value for said receiver, at said second selected time.
 14. Apparatus for estimating location coordinates of a signal receiver from signals received from two satellites, the apparatus comprising:a location determination unit that: receives a timed signal from each of first and second satellites at a receiver at approximately a first selected time; determines at least one of pseudorange value, Doppler shift value, a first selected combination of Doppler shift value and first time derivative of Doppler shift value, and a second selected combination of Doppler shift value and first and second time derivatives of Doppler shift value, with these values being collectively denoted location signal measurement values or LSM values, for the timed signal received from the first satellite; anddetermines an LSM value for the timed signal received from the second satellite; and a computer programmed to use the LSM values from the satellites to estimate first and second location coordinates values for the receiver, at a second selected time, determined with reference to the first selected time.
 15. The apparatus of claim 14, wherein:said location determination unit determines a selected LSM value for said timed signal for a selected one of said satellites at approximately said first selected time; and said computer is further programmed to estimate at least one of a receiver time offset value and a third receiver location coordinate value from the selected LSM value and said estimated first and second location coordinate values for said receiver.
 16. The apparatus of claim 14, wherein:said location determination unit determines a second LSM value for said timed signal for each of said first and second satellites at approximately said first selected time, where at least one of said first LSM values and the second LSM values for said two satellites is a pseudorange value; and said computer is further programmed to estimate a receiver time offset value and a third location coordinate value for said receiver from the second LSM values and said estimated first and second location coordinate values for the receiver. 